Patrones puntuales

library(RSocrata)
library(tidyverse)
library(lubridate)
library(tidygeocoder)
library(rgeoboundaries)
library(sf)
library(tmap)
library(spatstat)
source("Funciones_propias_ppp.R", encoding = "UTF-8")

Contextualización

La accidentalidad es uno de los problemas más comunes en el planeta provocando perdidas que pueden ir desde materiales hasta la vida misma, por lo que es de particular interés tratar de comprender como se comporta este fenómeno para evitar ser victimas de él.

La problemática en cuestión tiene la complicación de ser algo de naturaleza aleatoria, sin embargo esto no quiere decir que no exista alguna forma de controlarla pues la experiencia ha mostrado que las diferentes ciudades tienen lugares con mayor concentración de accidentes que otras ubicaciones dentro de la misma ciudad.

Debido a que se está trabajando la ocurrencia de accidentes en una ciudad de interés, la distribución Poisson y los procesos puntuales Poisson resultan ser una herramienta adecuada para tratar de explicar el fenómeno.

El estudio de accidentalidad se hará en los años 2019, 2020 y 2021 con el propósito de ver como se comporta la accidentalidad prepandemia, pandemia y postpandemia.

Obtención de la base de datos

Una vez presentado el fenómeno que se trabajará, se selecciona la ciudad de Barranquilla para realizar el estudio.

Los datos se extrajeron de la página web GOV.CO la cual tiene diferentes bases da datos de Colombia abiertas al público. ( Accidentalidad en Barranquilla). La base de datos contiene una gran cantidad de variables, sin embargo solo se escogen aquellas de interés para el patrón puntual. La estructura de los datos es presentada a continuación:

accidentes <- read.socrata("https://www.datos.gov.co/resource/yb9r-2dsi.json") %>%
    mutate(fecha = ymd(fecha_accidente)) %>%
    select(12, 6:8) %>%
    filter(year(fecha) > 2018 & year(fecha) < 2022)
knitr::kable(head(accidentes, 10),
             col.names = c("Fecha", "Gravedad del accidentes",
                           "Clase de accidentes", "Sitio accidente"))

Fecha	Gravedad del accidentes	Clase de accidentes	Sitio accidente
2019-01-01	Solo daños	Choque	CALLE 82 CRA 71
2019-01-01	Con heridos	Choque	CL 30 CR 27
2019-01-01	Con heridos	Choque	CL 17 CR 30
2019-01-01	Solo daños	Choque	CRA 31 CALLE 68C
2019-01-01	Con heridos	Atropello	CL 99D CR 8E
2019-01-01	Solo daños	Choque	VIA 40 CR 67B
2019-01-02	Solo daños	Choque	AV CIRCUNVALAR FRENTE AL NO 70A 12
2019-01-02	Con muertos	Choque	AV CIRCUNVALAR CR 13
2019-01-02	Solo daños	Choque	AV CIRCUNVALAR 49 39
2019-01-02	Solo daños	Choque	CR 53 68 216

Geocodificación de las direcciones

Se hace necesario realizarle unos ajustes a la base de datos considerada puesto que está no incluye las ubicaciones exactas de los accidentes ya sea en longitud - latitud o proyección UTM, para dicho propósito se usa la función geo() del paquete tidygeocoder para las direcciones de los accidentes en coordenadas de longitud - latitud.

direcciones <- paste(accidentes$sitio_exacto_accidente, ", Barranquilla, Colombia", sep = "")
localizaciones <- geo(address = direcciones, method = "arcgis")
write.csv(x = localizaciones, file = "accidentes.csv", row.names = F)

Área de estudio y base de datos final

bqlla <- geoboundaries(country = "COLOMBIA", adm_lvl = 2) %>%
    filter(shapeName == "DISTRITO ESPECIAL, INDUSTRIAL Y PORTUARIO DE BARR*")%>%
    st_transform(crs = 3857)

coordenadas <- read.csv("accidentes.csv")
coordenadas <- cbind(accidentes, coordenadas[, 2:3]) %>%
                st_as_sf(coords = c('long', 'lat')) %>%
                st_set_crs(value = 4326) %>%
                st_transform(crs = 3857) %>%
                st_intersection(bqlla)

accidentes <- list()
for(i in 2019:2021){
  accidentes[[as.character(i)]] <- filter(coordenadas, year(fecha) == i & month(fecha) > 3)
}

Gráfico de las localizaciones

Accidentes 2019

tmap_mode('view')

tm_shape(bqlla)+
  tm_polygons(alpha = 0.3, border.alpha = 0.7)+
  tm_shape(shp = accidentes[['2019']])+
  tm_dots(size = 0.01)

Accidentes 2020

tmap_mode('view')

tm_shape(bqlla)+
  tm_polygons(alpha = 0.3, border.alpha = 0.7)+
  tm_shape(shp = accidentes[['2020']])+
  tm_dots(size = 0.01)

Accidentes 2021

tmap_mode('view')

tm_shape(bqlla)+
  tm_polygons(alpha = 0.3, border.alpha = 0.7)+
  tm_shape(shp = accidentes[['2021']])+
  tm_dots(size = 0.01)

Pruebas de homogeneidad en la intensidad

Uno de los parámetros críticos al momento de modelar un parón puntual es la intensidad la cual puede ser constante (homogenea) o variable (inhomogenea), por lo tanto es importante iniciar el análisis con una prueba de homogeneidad de la intensidad.

# Guardando datos por anio
datos_2019 <- ppp(x = st_coordinates(accidentes[[1]])[, 1],
             y = st_coordinates(accidentes[[1]])[, 2],
             window = as.owin(W = bqlla))

datos_2020 <- ppp(x = st_coordinates(accidentes[[2]])[, 1],
             y = st_coordinates(accidentes[[2]])[, 2],
             window = as.owin(W = bqlla))

datos_2021 <- ppp(x = st_coordinates(accidentes[[3]])[, 1],
             y = st_coordinates(accidentes[[3]])[, 2],
             window = as.owin(W = bqlla))

Argumento gráfico

Conteos 2019

qc_2019 <- quadratcount(datos_2019, nx = 5, ny = 5)
plot(qc_2019)

Conteos 2020

qc_2020 <- quadratcount(datos_2020, nx = 5, ny = 5)
plot(qc_2020)

Conteos 2021

qc_2021 <- quadratcount(datos_2021, nx = 5, ny = 5)
plot(qc_2021)

Pruebas \(\chi^2\)

Accidentes en 2019

quadrat.test(qc_2019)

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 2857.1, df = 18, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 19 tiles (irregular windows)

Accidentes en 2020

quadrat.test(qc_2020)

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 1324.7, df = 18, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 19 tiles (irregular windows)

Accidentes en 2021

quadrat.test(qc_2021)

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 2094.2, df = 18, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 19 tiles (irregular windows)

Mapas de probabilidad

Año 2019

graph_ppp(datos_2019, 350, "Año 2019")

Año 2020

graph_ppp(datos_2020, 350, "Año 2020")

Año 2021

graph_ppp(datos_2021, 350, "Año 2021")